Self-consistent GW: an all-electron implementation with localized basis functions 
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This paper describes an all-electron implementation of the self-consistent GW (sc-GW) approach 
- i.e. based on the solution of the Dyson equation - in an all-electron numeric atom-centered or- 
bital (NAO) basis set. We cast Hedin's equations into a matrix form that is suitable for numerical 
calculations by means of i) the resolution of identity technique to handle 4-center integrals; and ii) a 
basis representation for the imaginary-frequency dependence of dynamical operators. In contrast to 
perturbative Go Wo, sc-GW provides a consistent framework for ground- and excited-state properties 
and facilitates an unbiased assessment of the GW approximation. For excited-states, we benchmark 
sc-GW for five molecules relevant for organic photovoltaic applications: thiophene, benzothiazole, 
1,2,5-thiadiazole, naphthalene, and tetrathiafulvalene. At self-consistency, the quasi-particle ener- 
gies are found to be in good agreement with experiment and, on average, more accurate than Go Wo 
based on Hartree-Fock (HF) or density-functional theory with the Perdew-Burke-Ernzerhof (PBE) 
exchange-correlation functional. Based on the Galitskii-Migdal total energy, structural properties 
are investigated for a set of diatomic molecules. For binding energies, bond lengths, and vibrational 
frequencies sc-GW and GoWo achieve a comparable performance, which is, however, not as good 
as that of exact-exchange plus correlation in the random-phase approximation (EX+cRPA) and 
its advancement to renormalized second-order perturbation theory (rPT2). Finally, the improved 
description of dipole moments for a small set of diatomic molecules demonstrates the quality of the 
sc-GW ground state density. 



Many-body perturbation theory (MBPT)± in the GW 
approach for the electron self-energy^— provides a natu- 
ral framework for an ah initio, parameter-free description 
of photo-ionization processes and charged excitations^ 
In recent years, the GW approach has become a popular 
method for the computation of band gaps and charged 
excitation energies for extended^ and finite systems^. 
In numerical implementations, following Hybertsen and 
Louie^ it is standard practice to treat the GW self- 
energy as a single-shot perturbation (Go Wo) acting on a 
Kohn-Sham (KS) or Hartrcc-Fock (HF) reference system. 
Thus, excitation energies are evaluated from first-order 
Feynman-Dyson perturbation theory as corrections to a 
set of single-particle eigenvalues. 

The popularity of the Go Wo approximation stems from 
the substantial reduction in the complexity of Hedin's 
equations at first-order perturbation theory: the KS or 
HF eigenstates from a self-consistent field calculation can 
be used as basis functions and provide a convenient rep- 
resentation in which the non-interacting Green function 
is diagonal. In this basis, only diagonal matrix elements 
of the self-energy £ arc needed to evaluate quasi-particle 
corrections at first-order. Thus, Go Wo grants a consider- 
able simplification of the linear algebra operations which 
is decisive for applying the theory to large molecules and 
solids. 

Although numerically more efficient than a non- 
perturbative approach, Go Wo suffers from several un- 
desirable shortcomings such as the dependence on the 
starting point^i^— the violation of conservation laws 
for momentum, total energy and particle number^ - — 



and - most importantly - the limited access to ground- 
state properties, which are kept unchanged from the pre- 
liminary density functional theory (DFT) or HF calcula- 
tions. 

It is known that the self-consistent GW approach (sc- 
GW) - in which both the Green function G and the 
screened Coulomb interaction W are iterated to self- 
consistency - ameliorates most of the pathologies of per- 
turbative GoWoJ^ A particularly appealing feature of 
the sc-GW method consists in the possibility of treating 
ground- and excited-states at the same level of theory. 
This property arises by virtue of the non-perturbative 
nature of the sc-GW approach, whereby the Green func- 
tion is updated and encompasses many-body effects in- 
troduced by the self-energy. In contrast, in perturba- 
tive theories (which generally do not introduce updates 
in the Green function) the electronic structure coincides 
with that of the corresponding starting point. Therefore, 
density and total energy - and derived quantities such 
as dipole moments, bond lengths, and binding energies - 
become accessible at self-consistency and reveal the qual- 
ity of the GW ground-state. Finally, at self-consistency 
excited- and ground-state properties are independent of 
the starting point, at least for closed shell systems^ and 
provide an unbiased assessment of the GW approach. 

A previous study on the homogeneous electron gas 
(HEG) reported a deterioration of the sc-GW spectral 
properties, as compared to GqWq*^ This has been at- 
tributed to a poor description of the satellite peaks 
at self-consistencyJ^ For extended systems, the perfor- 
mance of sc-GW remains controversial due to the scarce 



number of calculations for real solidsJ^r— Part of this 
controversy can be traced back to basis set problems 
in early all-electron calculations^ and to the large in- 
fluence that pseudo-potentials may have on GW band 
gaps^ More recently, sc-GW calculations for atoms^i 
and molecules^^ have shown improvements in the de- 
scription of the first ionization energies and for transport 
properties^ of finite systems. 

The price to pay in sc-GW is the demanding itera- 
tive procedure. The higher complexity of sc-GW arises 
for the following reasons, i) The Green function obtained 
from the solution of the Dyson equation is in general non- 
diagonal. This considerably increases the computational 
cost of the evaluation of the dielectric matrix, ii) The 
non-diagonal matrix elements of E are needed to solve 
the Dyson equation, iii) Fourier transforms of dynamical 
quantities are needed that introduce their own computa- 
tional difficulties. 

In the first part of this paper, we present an all-electron 
implementation of the sc-GW method in the localized 
basis-set code FHI-aims2£ and propose a recipe to ef- 
ficiently address points i) —iii) . An optimized set of lo- 
calized basis functions was used to represent the Green 
function and the self-energy operator. Non-local two- 
particle operators, such as the screened Coulomb inter- 
action, were computed by means of the resolution of 
the identity technique^— (also known as density fitting 
method) in a general framework previously introduced 
by some of us^I Finally, an auxiliary basis of Lorcntzian 
functions was introduced for an efficient analytical evalu- 
ation of Fourier transforms between imaginary time and 
frequency 

The second part of the paper, focuses on the assess- 
ment of ground- and excited-state properties as obtained 
from sc-GW for molecules. The quality of the sc-GW 
ground state was investigated by computing binding en- 
ergies, bond lengths, vibrational frequencies, densities, 
and dipole moments for a small set of hetero- and homo- 
atomic dimers. The full valence excitation spectra were 
evaluated for a set of molecules relevant for organic pho- 
tovoltaic applications (thiophene, benzoithiazole, 1,2,5- 
thiadiazolc, naphthalene, and tetrathiafulvalene) . From 
this study we conclude that sc-GW systematically im- 
proves the spectral properties of finite systems over the 
entire excitation spectrum (that is, not only for the first 
ionization energy) as compared to standard perturba- 
tive Go Wo calculations based on semi- local DFT and 
HF. Nonetheless, for certain starting points - exempli- 
fied by the PBEO hybrid functional - Go Wo slightly out- 
performs sc-GW, providing ionization energies in better 
agreement with experimental reference data, as also pre- 
viously demonstrated for benzene and the azabenzenes 
in Ref. Qjj. For structural properties the sc-GW method 
yields a less satisfactory agreement with experiment. For 
dimers, bond lengths and binding energies are slightly 
underestimated, and in this case there is no substantial 
improvement over perturbative approaches such as Go Wo 
or the random-phase approximation (RPA). Finally, self- 



consistency gives an accurate description of the electron 
density as manifested by the accurate dipole moments of 
diatomic molecules. These results suggest that sc-GW is 
a promising method for charge transfer compounds and 
interfaces. However, our study also indicates the impor- 
tance of including higher order exchange and correlation 
diagrams beyond GW to accurately describe the struc- 
tural properties of molecules. 

The paper is organized as follows: Section U gives a 
brief introduction to the GW approximation recalling 
the basic equations needed for the computation of the 
Green function G and the self-energy E. An optimal 
representation of Hedin's equations in terms of localized 
basis functions and of the resolution of the indentity is 
presented in Sec. [TIJ In Scc. lIIII we present the scheme em- 
ployed in the computation of the Fourier integrals of the 
Green function and other dynamical quantities. We re- 
port in Scc. lIVI an assessment of sc-GW for the excitation 
spectra of molecules and, in Sec. [V] for the ground-state 
properties of diatomic molecules. Our conclusions and 
final remarks are reported in Sec. IVI1 



I. THEORETICAL FRAMEWORK 

In MBPT the complexity of the many-body problem 
is recast into the calculation of the single-particle Green 
function. Knowledge of the Green function grants im- 
mediate access to the (charged) single-particle excitation 
energies of the system, to the total energy and, more 
generally, expectation values of any single-particle op- 
erator. Green function theory is well documented in the 
literature^ and we recall here only the basic equations rel- 
evant for the GW approach, adhering to Hartree atomic 
units Ti = m e = e 2 = 1. 

For a system of non-interacting electrons described 
through a time-independent Hamiltonian, the Green 
function can be written explicitly in terms of the single- 
particle eigenstates -0™ ( r ) an d eigenvalues e a n : 



G ff (r,r', W ) = ^ 



<(r)<*(r') 



w - « - M) - ir l sgn(fi - e£) 



(1) 



where n and a refer to orbital and spin quantum num- 
bers, respectively. \x is the electron chemical potential, 
and r\ a positive infinitesimal. In practice, V'nl 1 *) an d e^ 
are generally obtained from the self-consistent solution 
of the Hartrcc-Fock or (generalized) KS equations. 

For interacting electrons, the Green function has to be 
evaluated by solving the Dyson equation: 

G ff (r,r' ) w) = GJ(r,r , ,w) + f dr 1 dr 2 G^(v,r 1 ,u;)x 



x [E CT (n, r 2 , u) + Av5(n)(J(ri - r 2 ) 
-< c (r 1; r 2 )]G*(r 2 ,r» . 



(2) 



Here Aw^ is the change in the Hartree potential ac- 
counting for density differences between Gq and G 7 , and 



v° c is the exchange-correlation part of the single-particle 
Hamiltonian corresponding to the non-interacting Green 
function Go- For example, if Go is the Hartree-Fock 
Green function, then v° c corresponds to the non-local 
exchange operator E x . Alternatively, for a KS Green 
function, v° c is the local exchange-correlation potential. 
The electron self-energy E CT encompasses all many-body 
exchange-correlation effects and therefore its practical 
evaluation requires approximations. Following Hedin 2 -^, 
the self-energy can be expanded in a pcrturbative series 
of the screened Coulomb interaction W, with the first- 
order term given by the GW approximation: 



E CT (r, r', r) = iG a {r, r' , t)W{t, r',r) 



(3) 



By virtue of the time translation invariance it suffices to 
express E in terms of time differences (r = t — if). The 
screened interaction W(r, r', u>) is in turn defined through 
another Dyson equation: 

W(r,r',uj) = v(r,r')+ 

dr 1 dr 2 v(r 1 r 1 )x{ri 1 r 2l uj)W(r2,r',uj) . (4) 

Here, v is the bare Coulomb interaction l/|r — r'| and 
X the irreducible polarizability, which in GW is approx- 
imated by the product of two Green functions: 



X 



(r, r', t) = -iJ2 G CT (r, r', r)G CT (r', r, -r) 



(5) 



The self-consistent nature of Eqs. [2][5] arises from the 
interdependence of the self-energy and the Green func- 
tion. In the Go Wo approach, the self-energy is evaluated 
non-self-consistently, and the Dyson equation is solved 
approximately in a parturbative fashion. The quasi- 
particle excitation energies are then obtained from the 
quasi-particlc equation: 



f QP = e c 

n.cr r 



Re(C|S CT ( £ ^) 



Wu) 



(0) 



In this work, Eqs.[2][5]are solved fully self-consistently. 
In practice, an iterative procedure requires the following 
steps: 

1. Construction of an initial non- interacting Green 
function Go from a preliminary SCF calculation 
through Eq. [TJ 

2. Evaluation of the polarizability x from Eq. [5] and 
Fourier transformation of x to the frequency do- 



3. Calculation of the screened Coulomb interaction W 
from Eq. |4] and Fourier transformation of W to the 
time domain. 

4. Evaluation of the self-energy E from the Eq. [3] and 
Fourier transformation of E to the frequency do- 



5. Update of the Green function from the Dyson equa- 
tion (Eq. [2]) and Fourier transformation of G to the 
time domain. 

6. Mixing of the Green function to accelerate the con- 
vergence of the self-consistent loop. 

7. Iteration of steps 2.-5. until a convergence criterion 
is satisfied. 

In a numerical implementation, a choice for the basis 
set expansion of the quantities in Eqs.[2][5]has to be made. 
Our choice will be discussed in the next Section. 



II. SELF-CONSISTENT GW WITHIN A 
LOCALIZED BASIS 

Previous implementations of sc-GW were based on 
Gaussians or Slater orbitals , 18 ' 25 full potential linear aug- 
mented plane- waves , 20 ' 21 real-space grids,— and numeric 
atom-centered orbitals (NAO)^ In the present work, 
the Green function G, the self-energy E, and all single- 
particle operators, are expanded in a numeric atom- 
centered orbital basis {<fi(r)}, with basis functions of the 
form: 



<Pi(*) 



Ui(r) 



Y lm (n) 



(7) 



where Uj(r) are numerically tabulated radial functions 
and Y[ rn (Q) spherical harmonics. For numerical conve- 
nience, we work with real-valued basis functions by re- 
quiring - without loss of generality - that Yj m (57) denotes 
the real part (for m — 0, . . . , /) and the imaginary part 
(for m = — /,...,— 1) of complex spherical harmonics. In 
FHI-aims the choice of the radial functions u;(r) is not 
limited. In this work we will show results for numerically 
tabulated Gaussian orbital basis sets and the Tier hier- 
archy of FHI-aims for NAOs^i We refer to Ref. H3 for 
details on the construction and optimization as well as 
the properties of the NAO basis sets in FHI-aims. 

In terms of the <fi basis functions, the Green function 
can be expanded as: 



A',, 



G ff (r, 



£ p i (r)*r. 1 G$M«r m Vm(r') , (8) 

ijlm 



where s^ = J dripi(r)ipj(r) is the overlap matrix tak- 
ing into account the non-orthonormality of the basis set 
and iVbasis is the total number of basis functions. In the 
following, sums over latin indixes i,j,l,m are implicitly 
assumed to run from 1 to iVbasis ; whereas sums over n run 
over the total number of states. The coefficients GfAiui) 
of the expansion are given by: 



Gf j (iu)= J drdrVi(r)G CT (r,r',w)^(r') 



(9) 



The representation in Eq. [9] can be easily applied to the 
non-interacting Green function in Eq. [TJ yielding 



W=EE 



S H C ln C mn s mj 
i-U ~ (< ~ (J.) 



(10) 



where we introduced the expansion of the HF/KS orbitals 
in the NAO basis Vn( r ) = Sz c to<M r ) an d the Green 
function was continued to the imaginary frequency axis. 
The matrix representation in Eq. [§] is also adopted for 
the Hartree potential ug, the self-energy S cr , and the 
exchange-correlation potential v° c . 

To rewrite Hcdin's equations in a matrix form suitable 
for a numerical implementation, we need to introduce 
a matrix representation for two-particle operators. The 
expansion of two-particle operators in a numerical basis, 
cannot be handled efficiently through Eq. [5] due to the 
appearance of the 4-orbital 2-electron Coulomb integrals 
of the form: 



(ij\M) 



tpi(r)ip j (r)tp k (r')<p l (r') 



drdr' 



r - r' 



(11) 



The computation of the Coulomb repulsion integrals 
in Eq. [TJJ is a problem extensively discussed in the 
literature^— and efficient techniques have been devel- 
oped over the years to make this calculation affordable. 
Numerically, the difficulty arises from the large number 
of NAO pairs and from the memory requirements of stor- 
ing the 4- index matrix (ij\kl). In the NAO framework, 
the integrals in Eq.[TT]are often evaluated by introducing 
an auxiliary basis set {P M (r)}, with basis functions -P^(r) 
defined to span the Hilbcrt space of NAO pairs 



ViMVjM 



m=i 



C£P„(r) 



(12) 



where C£ are the coefficients of the expansion. Due to 
the high linear dependence of the NAO products, the 
number of product basis functions iV aux is much smaller 
than the number of NAO pairs 0(N^ asis ), making the 
numerical evaluation of Eq . 1 1 1 1 affordable . This techique, 
known as the resolution of the identity (RI) - or also 
density-fitting technique - was implemented in the FHI- 
aims code and we refer to Ref. [3l| for a detailed account of 
the variational approach employed in the determination 
of the RI coefficients Cjt- and for a review of the overall 
accuracy of the RI approach for correlated calculations. 
In short, we used the "RI-V" variant of the RI scheme, 
in which the expansion coefficients are given by: 



C£ = X>» 



V, 



1/fJ, 



(13) 



where (ij\v) = Jdr< y 9 l (r)^ i (r)P /i (r')/|r-r / | and V vii de- 
notes matrix elements of the Coulomb matrix in the aux- 
iliary basis, i.e., V^ — J drdr'P^(r)P^(r')/\r — r'|. For 
numerical efficiency, it is convenient to work with the 
generalized coefficients: 



M£ 






/2 



(14) 



Following Ref. |3l|, one can write the RI version of 
the Dyson equation for the screened Coulomb interac- 
tion (Eq. gj as: 

WVM = [v^wiiu)]^ = [i - n(* w )]^ , (is) 

where we defined [II(iu;)] = [x(«w)i>] . In contrast 
to the Rl-based implementation of the Go Wo method^!, 
the operator II has to be expressed as an explicit func- 
tional of G. Moreover, all non-diagonal matrix elements 
in the Green function have to be included. These two 
criteria are satisfied by evaluating II in terms of the Mf- 
coefficients, as: 

P(»t)] m „ = ^Y,H Ki^mGl^Glmi-ir). (16) 

<7 ijlm 



Here we defined 

Gij(ir) = J2 s a lG im{iT)s;^ 



(17) 



The quadruple sum in Eq. [16] may be reduced to dou- 
ble sums - with a considerable reduction of computa- 
tional cost - by introducing the intermediate quantity 
Af- CT (ir) = Y^,i ^ifiij^ T )- l n terms of these coefficients 
Eq. [TBI becomes: 



mr)U = -i^2^2A^ a (ir)A^ a (-ir) 



(18) 



The self-energy can be evaluated in terms of Eq. [15] 
providing the following matrix representation of Eq. [3] 



E W = ^EE M « M ^m( ir ) W ^ ir ) . (19) 

lm fiv 



By introducing the auxiliary quantity B^ m (ir) = 

Y^u-^jmWfi"(^ T )> the self-energy can again be cast into 
a double-sum form: 



XUir) 



2tt 



EE^ )ff (*r)^ m (iT) 



(20) 



The correlation (exchange) contribution to the self- 
energy can be derived straightforwardly from Eq. [19] by 
substituting W with W^ = W ^ - 6^ (W*„ = 6^ v ). 
The Hartree potential v^ is also evaluated as an explicit 
functional of G as: 

^ = J2J2K M L^m(ir = 0-) . (21) 

lm ^ 

Finally, the matrix representation of the Dyson equation 
for the Green function completes the set of Hcdin's equa- 
tions: 



G„-M= GoM^-xtm + c 



A« H 



-l 

ij 



(22) 
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FIG. 1. (Color online) Panel (a): values of A - defined in 
Eq. [24] - as a function of the number of iterations of the sc- 
GW loop for H2, H2O and CeHe in their equilibrium geometry 
in a Tier 2 basis set. A linear mixing parameter a — 0.2 was 
used for all molecules. Ath indicates the default value of the 
convergence threshold. Panel (b): values of A for H2O as a 
function of the number of iterations for different values of a. 
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FIG. 2. (Color online) Total time (in seconds, on a single 
CPU) per iteration of the sc-GW loop, for linear hydrogen 
chains of different lengths. The total time required for the 
evaluation of the self-energy in GqWq is included for compar- 



Here Aug is the difference between the Hartree potential 

of the interacting and the non- interacting Green function. 

To facilitate the convergence of the sc-GW loop, the 

input Green function G of the (N + l)-th iteration 

is obtained from a linear mixing scheme: 



G ij (« T ) 



■nN 



a&iAiT) + (I - aj&ij (it) 



(23) 



where G denotes the Green function obtained from the 
iV-th solution of the Dyson equation, and a is the mixing 
parameter. As illustrated in panel (b) of Fig. [TJ we find 
that a = 0.2 is typically a good choice. The convergence 
of the self-consistent loop is monitored looking at the 
average deviation of the Green function at each iteration 
A, defined as: 



1 



basis 



E 

i,3 



|Gy (it = 0") - G£ V = 0- 



(24) 



The sc-GW loop is considered converged when A drops 
below a chosen threshold A t h- Default settings used in 
most calculations are A t h = 10~ 5 , which suffices to con- 
verge both total and quasi-particle energies. The con- 
vergence of sc-GW is illustrated in Fig. [TJ where A is 
reported as a function of the number of iterations for H2 , 
H 2 and C 6 H 6 . 

Equations I15H22I constitute a matrix representation of 
Hcdin's equations in the GW approach (Eqs. [5][5]). We 
emphasize again that in Eqs. [15l[22l i) all electrons are 
treated on the same quantum mechanical level, i.e. fully 
self-consistcntly; ii) no model screening was used in the 
calculation of W; iii) all non-diagonal matrix elements of 
G and S are correctly accounted for. 

The evaluation of Eqs. \W\ and [19] is the most compu- 
tationally demanding operation of our implementation. 
The scaling of the computational time as a function of the 



basis set size equals that of Go Wo calculations but with 
a larger prefactor. To illustrate this aspect, we report 
in Fig. [2] the total computational time spent for a single 
iteration of Eqs. [T5l[22l as function of the length of a lin- 
ear hydrogen chain in a minimal basis set (i.e., with one 
NAO per atom). As compared to conventional Go Wo im- 
plementations, the additional computational cost arises 
from the necessity of accounting for non-diagonal matrix 
elements in the calculation of G and S. 



The only approximation introduced up to this point 
is the resolution of the identity for the expansion of the 
product of NAO pairs (Eq.H2]). As discussed in Ref. HH, 
the accuracy of the RI can be monitored systematically 
by means of two control parameters: £ or th and £svd- 
Sorth sets the accuracy threshold for the Gram-Schmidt 
orthonormalization employed for the reduction of the lin- 
ear dependence of on-site (i.e. on the same atom) product 
basis functions P^. In practice, by chosing smaller values 
of £ rth one may increase the number of product basis 
functions used in the expansion in Eq. 1121 Similarly, the 
parameter £svd controls the singular value decomposi- 
tion (SVD) for the orthonormalization of product basis 
functions on different atoms. A more detail description 
of the effects of these parameters can be found in Ref. [3X|. 
To show the effect of the RI scheme on the self-consistent 
Green function, we report in Fig. [3] the sc-GW total en- 
ergy - evaluated from Eq. [3D] introduced in Sec. [VI— of 
the water molecule as a function of e rth (left panel), and 
£svd (right panel). For a wide range of values of the 
control parameters £ rth and £svd, the changes of the 
total energy are of the order of 10" 4 cV or less. In all 
following calculations we therefore used e rth = 10~ 5 and 
esvd = 10" 5 . 
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FIG. 3. sc-GW total energy of H2O as a function of the con- 
vergence parameters e or th (left panel) and esvd (right panel) , 
evaluated with a Tier 2 basis set. The number of product 
basis functions corresponding to each value of Eorth is also 
reported. 
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FIG. 4. (Color online) Comparison of the analytic structure of 
the real (in blue, above) and imaginary part (in orange, below) 
of f n (ioj) for different values of b n , and a matrix element of 
the Green function (black). 



III. DISCRETIZATION OF THE FOURIER 
INTEGRALS 

In this implementation, we solve Eqs. I15H22I in imagi- 
nary time and frequency, taking advantage of the reduced 
number of frequency points required to describe G° (iui) 
and other dynamical quantities, as compared to real fre- 
quency implementations. In a mixed time-frequency for- 
malism convolutions on the frequency axis can be ex- 
pressed, by virtue of the convolution theorem, as prod- 
ucts on the time axis after a Fourier transform. Due to 
the slow decay of G a {ioj) at large frequencies, the com- 
putation of Fourier transforms may require extended and 
dense frequency grids. We obviate this problem by intro- 
ducing a basis for the frequency /time dependence of all 
dynamical quantities. This permits an analytic evalua- 
tion of the Fourier integrals - as one can choose basis 
functions with a Fourier transform known analytically - 
and substantially reduces the number of frequency points 
needed to converge the calculation. 

Following the approach introduced in Ref. |4Gj, we ex- 
pand the Green function in a set of Lorentzian functions 
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FIG. 5. (Color online) Mean absolute error (MAE) introduced 
by Fourier transforming the Green function of N2 for different 
numbers of poles and frequency points. 



of the form f n (iu) = l/(b n + iui), with Fourier transform 
f(ir) = l/(27r)e~ 6 " T . The parameters b n are fixed at the 
beginning of the calculation and are distributed logarith- 
mically in the energy range covered by the Kohn-Sham or 
Hartree-Fock eigenvalues e°- Although in principle other 
functions could be used, the functions f n {iui) constitute a 
natural choice for the expansion of the Green function, as 
the frequency dependence of f n {iui) closely resembles the 
analytic structure of G and captures the \/iui behaviour 
at large frequency. This is illustrated in Fig.EJ where the 
real and imaginary parts of f n (iu>) - with different values 
of b n - are compared to a Green function matrix element 
for the Ne atom. The Green function can be expanded 
in the basis of Lorentzian functions as: 



A',, 



G"M= £ K^H + ^/TH] 



(25) 



where iV po i 0S denotes the number of functions f n , and 

fn e (iu) the real (imaginary) part of f n (iuj). The 
real and imaginary part of the Green function have been 
treated separately to maintain a real-valued linear-least 
square problem, leading in turn to real- valued coefficients 
a n and /?„. Since the Fourier transform of the f n {iu)) is 
known, the coefficients a n and /?„ also determine the ex- 
pansion of the Green function in imaginary time. Expan- 
sions similar to Eq. l25l wcrc employed also for the Fourier 
transform of Xi W an d S. 

The imaginary time and frequency axes are then dis- 
cretized on exponentially spaced grids composed of N^ 
points in the range {0,w ma x}, and by 2N T + 1 points in 
the range {— T max , r max }. The grid points uj n and inte- 
gration weights w(u) n ) are defined as: 



U>k 



loo [e (fc - 1)h - l] , w(uj k ) = hujoe (k -^ h (26) 



and similarly for r^ and w(rk). The constant h is ob- 
tained by imposing the maximum frequency w max from 



the constraint a; n 



UJQ [( 



N^h 



l] and the parame- 



tcr ljo sets the initial spacing of the grid. Typical val- 
ues adopted in our calculations are w max = 5000 Ha, 
r max = 1000 Ha" 1 and uj q = r = 0.001. 

The error introduced by the Fourier transform can 
be quantified for functions known analytically on both 
the (imaginary) frequency and time axes such as, for in- 
stance, the non- interacting Green function given in Eq.[T] 
In Fig. [3J we report the mean absolute error (MAE) in 
the Fourier transform of the non-interacting Green func- 
tion Gq(ilu) of the nitrogen dimer N 2 , averaged over all 
matrix elements. The MAE drops exponentially when 
increasing the number of functions /„, and few tens of 
frequency points suffice to converge the Fourier integrals 
with an accuracy of the order of 10~ 8 . In our calculations 
we used N u = N T = 60 as default parameters. 



IV. SPECTRAL PROPERTIES OF MOLECULES 

We turn now to the spectral properties in sc-GW. At 
self-consistency, the excitation spectrum is given by the 
spectral function: 

A(oj) = —1/it / dr lim Im G(r, r',w) 

J r '^ r 

= -1/tt Tr [Im G(u)] , (27) 

where the Green function has to be evaluated on the 
real frequency axis. To evaluate Eq. 1271 we first obtain 
the real frequency self-energy by means of the analytic 
continuation based on a two-pole fitting schemed. In this 
approach, the matrix elements of the self-energy in the 
imaginary frequency domain (i.e., the Fourier transform 
of Eq. [20)) are fit by polynomials of the form: 



S(zw) 



2 

E 



b„ 



(28) 



Here the matrix element indices were suppressed for no- 
tational simplicity and the coefficients a n and b n are de- 
termined by means of a non- linear least-square fit, solved 
with a Lcvcnbcrg-Marquardt algorithm. By replacing iu> 
by uj in Eq. [55] the self-energy can then be evaluated 
on the real frequency axis. Once the real-frequency self- 
energy is obtained, the Dyson equation is solved directly 
in real frequency on a fine, equally spaced grid. The re- 
sulting Green function is used to determine the sc-GW 
spectral function A{uS). 

Previous works^i have indicated that the two-pole 
model presented in Eq. [2FJ reliably reproduce quasi- 
particle energies with an average relative error of the 
0.2% for solids. The parameter r\ in the denominator 
of Eq. [TJ necessary to avoid the numerical divergence of 
Go is set to r\ = 10~ 4 . This parameter contributes negli- 
gibly to the broadening of the spectral function and has 
no effect on the position of the quasi-particle peaks. 

As an example, we report the sc-GW spectral func- 
tions of H2O, NH3 and N2 in Fig. [S] calculated using ba- 
sis sets of increasing size. The sc-GW spectral function 



shows sharp (5-function-like peaks at the quasi-particlc 
energies. The absence of broadening in the quasi-particle 
peaks in Fig. [5] may be associated with a infinite lifetime 
of the corresponding quasi-particle states, as expected 
for states close to the Fermi energy. As discussed in 
Sec. IIV Al higher energy excitations may decay through 
the formation of electron-hole pairs, leading to a finite 
lifetime and thus to a more pronounced broadening of the 
quasi-particle peaks. In panels a), b) and c) in Fig. [6] we 
report the spectral function corresponding to the high- 
est occupied quasi-particle states evaluated with a Tier 
1, Tier 2 and Tier 3 basis; panels d), e) and f) show 
the peaks corresponding to the lowest unoccupied quasi- 
particle states. The GoW / o@HF and sc-GW ionization 
energies are reported in panels g), h) and i) of Fig. [5] as 
a function of the basis set size. The GqWq ionization en- 
ergy is calculated from the linearized quasi-particle equa- 
tion (Eq. [5]), whereas in sc-GW it is extracted from the 
highest (valence) peak of the spectra shown in panels a) , 
b) and c). 

For the quasi-particle energies corresponding to occu- 
pied states, the largest change is observed going from 
Tier 1 to Tier 2. For N2 for example, we observe a change 
in the HOMO of approximately 0.2 eV going from Tier 
1 (which consists of 14 NAO basis functions per atom) 
to Tier 2 (39 NAO per atom). A further increase of 
the basis set from Tier 2 to Tier 3 (55 NAO per atom) 
leads to changes of the order of 5 meV in the HOMO 
as illustrated in the right panels of Fig. [5] Lower lying 
quasi-particle peaks show a similar convergence behav- 
ior as the HOMO. H2O, and NH3 exhibit a qualitatively 
similar behavior. Surprisingly, for all systems considered 
here, sc-GW data converge faster with the basis set size 
than perturbative GqWq calculations. In the following, 
we will focus on closed shell molecules, which in many 
instances do not have a stable anionic state. Therefore, 
we will focus on the spectral function of occupied states 
only. 

To investigate the performance of the GW approxima- 
tion at self-consistency, we have performed sc-GW calcu- 
lations for a set of molecules relevant for organic photo- 
voltaic applications. We report in Fig. [7] the comparison 
between experimental^— and theoretical ionization en- 
ergies evaluated from sc-GW and GqWq based on the 
HF, PBE, and PBE0 starting points, for thiophene, ben- 
zothiazole, 1, 2, 5-thiadiazole, naphthalene, and tetrathi- 
afulvalene. For an unbiased assessment, it would be de- 
sirable to benchmark sc-GW against higher level theo- 
ries, since in experiment the distinction between vertical 
and adiabatic ionization energies is difficult and vibra- 
tional effects are always present. For naphthalene the 
coupled cluster singles doubles with perturbative triples 
(CCSD(T)) method, that is currently considered as the 
gold standard in quantum chemistry, gives a vertical ion- 
ization potential of 8.241 eV— , which sc-GW underes- 
timates (-7.48 eV). For this molecule, the difference be- 
tween the vertical and the adiabatic ionization potential 
is only 0.1 cV in CCSD(T). For thiophene, CCSD(T) 
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FIG. 6. (Color online) Spectral function calculated from Eq. 1271 using a sc-GW Green function for H2O, NH3 and N2 with Tier 
1 (black), Tier 2 (orange) and Tier 3 (blue, dashed) NAO basis sets. Panels a), b) and c) show the peaks corresponding to the 
first valence states, peaks relatives to conduction states are reported in panels d), e) and f). Finally, the panels g), h), and i) 
report the convergence of the HOMO level as a function of the number of numerical orbitals used in the basis set. 



calculations of the adiabatic ionization energy obtain 
8.888 eV— , in good agreement with experiment, whereas 
sc-GW yields 8.45 eV. Zero-point vibration effects are 
small and cancel with core-correlation and relativistic ef- 
fects. However, the authors of this study indicate that 
the geometry of the cation differs considerably from that 
of the molecule, but did not give values for the vertical 
ionization energy. It therefore remains an open ques- 
tion, by how much vertical and adiabatic ionization po- 
tentials differ for thiophenc. For benzothiazole, 1,2,5- 
thiadiazolc, and tctrathiafulvalene we were not able to 
find CCSD(T) calculations for the vertical ionization po- 
tential. 

Despite the tendency to underestimate the first ion- 
ization energy, for these systems sc-GW ionization ener- 
gies are in good agreement with experiment, and give 
a good overall description of the excitation spectrum: 
full self-consistency leads to an average error of 0.4 eV 
(with a maximum error of 1.2 eV) between the experi- 
mental and theoretical ionization energies, whereas HF- 
and PBE-based GqWq differs on average by 0.7 eV (with 
a maximum error of 1.5 eV for GoWo@HF, and 1.6 eV 
for G W @PBE). Interestingly, G W @PBE0 ionization 
energies arc close to the sc-GW ones. Moreover, the 



GoWo@PBEO spectrum is in slightly better agreement 
with experiments with an average deviation of 0.3 eV 
(with a maximum error of 1.2 eV) - as recently also re- 
ported for benzene and the azabenzenes in Ref. Il2l . 

For small molecules, the improvements of the spectral 
properties at self-consistency can partially be traced back 
to corrections of the over- or under-screening in Go Wo- In 
PBE based G W calculations, the small HOMO-LUMO 
gap induces an overestimation of the screening in the 
Coulomb interaction. This is the origin of a systematic 
error in the GoW / o@PBE quasi-particle energies. Simi- 
lar considerations are easily generalized to the HF start- 
ing point, where HOMO-LUMO gaps are generally too 
large due to the missing correlation energy. PBE0, on the 
other hand, often gives a good compromise. Therefore, 
the over- and under-screening is reduced in GoWo@PBEO 
and the resulting excitation spectrum is similar to sc- 
GW, where - because of the self-consistent calculation of 
W - this problem is mitigated. Based on these results, 
PBEO emerges as an optimal starting point for the per- 
turbative calculation of the spectral properties. It was 
argued that the screened Coulomb interaction obtained 
from sc-GW may also be underscreened due to the lack 
of electron-hole interactions - typically accounted for by 
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FIG. 7. (Color online) Comparison between theoretical and 
experimental vertical ionization energies of thiophene, ben- 
zothiazole, 1, 2, 5— thiadiazole, naphthalene, and tetrathiaful- 
valene. Experimental photoemission data are from Refs. 1421 - 
l4a . The molecular geometries were optimized with PBE in a 
Tier 2 basis set and are reported on the right. Go Wo ioniza- 
tion energies are obtained with a Tier 4 basis set, the sc-GW 
ones with a Tier 2 basis. 



vertex corrections^. This would in principle lead to a 
systematic overestimation of the quasi-particle energies, 
as reported in Ref. [50 for semiconductors. Figure [3 on 
the other hand, indicates a slight underestimation of the 
sc-GW quasi-particle energies, confirming the expecta- 
tion that effects of the electron-hole attraction on the 



screened Coulomb interaction are small in molecules with 
large HOMO-LUMO gaps. 

As alluded to in the introduction, previous sc-GW 
studies have reported conflicting conclusions on the accu- 
racy of the spectral properties jL12 r 20 i 25 Consequently, no 
consensus has so far been reached in this respect. sc-GW 
calculations for the homogeneous electron gas (HEG) in- 
dicated a deterioration of the spectra as compared to 
perturbative Go Wo based on the local-density approx- 
imation (LDA)J^ For the HEG, Holm and von Barth 
observed a transfer of spectral weight from the plas- 
mon satellite to the quasi-particle peak in self-consistent 
calculations^ This results in a weaker plasmon peak and 
a broader valence band, that worsens the agreement with 
photo-emission experiments for metallic sodium. 

The first self-consistent calculation for real systems 
- performed for potassium and silicon in the pseudo- 
potential approximation - confirmed the picture outlined 
by Holm and von Barth, indicating a deterioration of the 
band width and band gap at self-consistencyJ^ In a later 
work, Ku and Eguiluz attributed the origin of this fail- 
ure to the pseudo-potential approximation, emphasizing 
the importance of accounting for core- valence coupling in 
the determination of the screening^. However, several 
groupSi 23 ' 51 i 52 have questioned the convergence of these 
calculations with respect to the number of bands. Nev- 
ertheless, these earlier studies gave the impression that 
full self-consistency deteriorates the spectral properties 
compared to perturbative Go Wo, and that it is not rec- 
ommended to perform sc-GW calculations. However, in 
our opinion, the scarce numerical evidence for realistic 
systems is not enough to corroborate this notion. 

It was argued that the deterioration of spectra in sc- 
GW might arise due to the iterative construction of the 
polarizability \ as the product of two Green functions^. 
This would systematically weaken the incoherent part of 
the Green function, and reduce the intensity of the plas- 
mon satellites. For molecules however, this mechanism 
does not apply since quasi-particle peaks carry integer 
spectral weight, and no plasmon satellites are observed. 
For extended systems, this mechanism might effectively 
deteriorate the sc-GW spectral function, as for the homo- 
geneous electron gas. Nonetheless, more investigations 
are needed to provide a general and systematic assess- 
ment of sc-GW for real solids. 



A. Lifetimes of quasi-particle peaks 

To facilitate the comparison between GoWo@PBEO 
and sc-GW, we report in Fig. [8] the full sc-GW and 
GoWo@PBEO spectral function for thiophene and 1,2,5- 
thiadiazole. Figure |S] illustrates that even if the peak 
positions in sc-GW and GoWo@PBE0 are very similar, 
there are qualitative differences. In GoWo@PBE0 the 
lower lying peaks are considerably more broadened than 
the corresponding sc-GW peaks. 

We observe that quasi-particle peaks corresponding 
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FIG. 8. (Color online) Comparison between sc-GW and GoWjQPBEO spectral functions for thiophene and 1, 2, 5— thiadiazole 
evaluated with a Tier 2 basis set. Vertical dashed lines indicates experimental photoemission data from Refs. |42| and 143 . 
respectively. GoVKqQPBEO spectra are obtained with a Tier 4 basis set, the sc-GW spectra with a Tier 2 basis. 



to high-energy excitations are accompanied by a finite 
broadening. The broadening in turn, being inversely 
proportional to the lifetime of the corresponding quasi- 
particle state, yields important information on the dy- 
namics and damping of excitations. A careful quanti- 
tative evaluation of lifetimes would require to properly 
account for effects beyond the GW approximation (such 
as the coupling between electrons and vibrations and the 
satisfaction of selection rules in the decay process) and 
goes beyond the purpose of the present work. There- 
fore, the following discussion will be limited to explain 
why even the quasi-particle excitations of molecules have 
finite lifetimes and we will briefly characterize their start- 
ing point dependence. 

The physical origin of the lifetime is simple. The hole 
created in a lower valence (or core) state by the photoe- 
mission process can in principle decay to the HOMO, or 
to states energetically close to the Fermi energy. The en- 
ergy released in this process has to be converted into an 
internal excitation of the system, since isolated molecules 
cannot dissipate energy. If the energy released is larger 
than the HOMO-LUMO gap, a particle-hole pair can be 
created. This opens up a scattering or decay channel for 
the hole, which therefore acquires a lifetime. The energy 
threshold for electron-hole formation is then given by 
A = £gg M o - £g G ap, with ££g MO and E^ p the HOMO 
level and the HOMO-LUMO gap of the starting point, 
respectively. In other words, only quasi-particle states 
with an energy below A may decay, and acquire a finite 
broadening. This argument is general and does not only 
apply to GW. What is particular to GqWq is that the 
relevant gap for this process is determined by the start- 
ing point; in this case the DFT functional for the ground 
state. 



To illustrate this effect on the broadening of the quasi- 
particle peaks, we report in Fig. [S] the spectrum of ben- 
zene evaluated from Go Wo based on different starting 
points and at self-consistency. Values of A for the dif- 
ferent exchange-correlation functionals are reported as 
vertical dashed lines (in green). PBE has the smallest 
HOMO-LUMO gap (~ 5.2 eV) and we observe a notice- 
able peak broadening (i.e., short lifetime) at ~14 eV in 
the GoW)@PBEO spectrum. A systematic increase of the 
broadening is then observed the further a state lies be- 
low A. Adding exact exchange to the DFT functional 
and increasing its admixture opens the HOMO-LUMO 
gap. The onset of the finite lifetime subsequently moves 
to lower energies. 

In sc-GW, the broadening of the quasi-particle peaks 
is consistent with the HOMO-LUMO gap at the GW 
level, and the ambiguity of the starting point depen- 
dence is lifted. Consequently, for benzene only peaks be- 
low — 19 eV acquire a small finite broadening. Based on 
the sc-GW results, the large broadening observed in the 
GoWo@PBEO spectrum can be attributed to the small 
HOMO-LUMO gap of the underlying PBE0 calculation. 
The inclusion of a fraction of exact exchange partially 
ameliorates this problem, but not fully. This indicates 
that the calculation of lifetimes presents an additional 
situation in which - due to the severe dependence on the 
starting point - resorting to full self-consistency is essen- 
tial. 



B. Independence of the starting point 

Previous work showed that partially self-consistent 
approaches such as eigenvalue self-consistent GW— or 
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FIG. 9. (Color online) Spectral function of benzene evaluated 
from sc-GW and GoWo based on PBE, PBEO with differ- 
ent mixture of exact exchange (EX), and HF. Vertical dashed 
lines (in green) indicate the energy threshold A for the for- 
mation of an electron-hole pair, which depends explicitly on 
the HOMO and LUMO levels of the underlying DFT/HF cal- 
culation. The quasiparticle peaks acquire a finite broadening 
only for states below A. For peaks above A, the residual 
broadening stems from the parameter r\ = 10 -4 discussed in 
the text. 



quasi-particlc sclf-consistcncy ) 53 i 54 reduce the starting 
point dependence but they do not eliminate iti 12 ' 55 Only 
full self-consistency successfully removes any dependence 
on the starting point, as we discussed in Ref.[l7l This is a 
major advantage of the sc-GW scheme, as it allows a sys- 
tematic assessment of the GW approximation unbiased 
by spurious dependence on the input Green function. 

To illustrate the independence of the starting point 
of the sc-GW Green function, we report in Fig. [TO] the 
spectral function of the carbon monoxide molecule as a 
function of the number of iterations of the Dyson equa- 
tion initialized with HF, PBE, and PBEO. After just a 
few iterations of the sc-GW loop, the quasi-particle peaks 
in the spectral function arc located at the same energies 
demonstrating the independence of the starting point in 
sc-GW. 



V. GROUND-STATE PROPERTIES FROM THE 
GW APPROXIMATION 

A. Galitskii-Migdal total energy 

In MBPT, the total energy E to t may be regarded as 
an explicit functional of the Green function, i.e. .Etot = 
Etot[G]. The functional dependence of E tot [G] on G is 
not uniquely defined and different total energy function- 
als have been proposed over the years. Some example are 
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FIG. 10. (Color online) Spectral function of the CO dimer 
evaluated for the first five iterations of the Dyson equation 
from the HF, PBE, and hybrid PBEO starting points. The 
first iteration corresponds to the GqWq approximation. 



the Luttinger-Ward^ and the Klein functional^ - which 
are stationary at the self-consistent Green function^ - 
and the Galitskii-Migdal formula^. The total energy 
obtained from different functionals may differ in princi- 
ple, if evaluated with a given Green function. However, 
if the Green function is obtained self-consistently solv- 
ing the Dyson equation, all functionals yield the same, 
unique total energy. Since we are interested in sc-GW^ 
total energies, all total energy functionals arc equivalent 
and we use the Galitskii-Migdal formula because of its 
simplicity: 



Egm — —i 



y2pTr{[u> + h ]G°(Lj)}+E ion . (29) 



Here /io is the single-particle term of the many-body 
Hamiltonian, i.e. the sum of the external potential due 
to the nuclei and the kinetic energy operator, and E[ on 
accounts for the repulsive nuclear energy. As discussed 
in Ref. [32|, Eq. [29] can be computed directly in imagi- 
nary frequency taking advantage of the reduced size of 
the integration grid needed to describe G and E on the 
imaginary axis. In the present work, we cast Eq. 1291 into 
a more suitable form for numerical implementations^: 



E, 



GM 



;£g;.(t = o-)[2^ + 2 



v*f +vf i + ^ i J 



t],o 



*E 



^•ME^Me^ + i-U, (30) 

Z7T J J 



where t denotes the kinetic-energy operator, w H and 
v ext the Hartree and external potential, and E x and E c 
the exchange and correlation parts of the self-energy, 
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FIG. 11. (Color online) BSSE corrected binding energy of CO 
evaluated with Tier 1, 2, 3 and 4. Each curve is aligned at 
their respective Tier 4 value. All data are in eV. 



respectively. We used G 



l Gs- 



with 



f dripi(r)ipj(r). A derivation ol Eq. [30] is reported in 
Appendix fX] We emphasize that Eq. [3UJ is exact if eval- 
uated with the exact self-energy and Green function. In 
Eq. 1301 contributions arising from time- independent oper- 
ators can be easily evaluated by simple matrix products. 
The correlation energy on the other hand requires a fre- 
quency integration which can be evaluated directly on 
the imaginary frequency axis without resorting to the 
analytic continuation 3 - 2 -. We evaluated Eq. [30] in the 
GW approximation. In sc-GW, both E and G are self- 
consistent solution of the Dyson equation, whereas for 
GqWq the self-energy is evaluated only once and G is 
the non-interacting Green function of the DFT/HF cal- 
culation. The latter procedure corresponds to a first- 
order perturbative correction of the DFT/HF total en- 
ergy, with a perturbing potential given by [E(w) — v xc ]. 



B. Structural parameters of diatomic molecules 

Total energy differences are more important than abso- 
lute total energies, as they give information on structural 
properties and on the strength of chemical bonds. Here, 
we restrict the discussion to the ground-state properties 
of dimers at their equilibrium geometry. The reader is re- 
ferred to Ref. [6l| for an assessment of the sc-GW method 
in the dissociation limit. 

In the following, we provide an assessment of the sc- 
GW method for bond lengths, binding energies, and vi- 
brational frequencies based on the calculations of the po- 
tential energy curve of LiH, LiF, HF, CO, H2, and N2. To 
illustrate the convergence with the basis set, we report 
in Fig. [TT] the binding energy of the carbon monoxide 
dimer evaluated with sc-GW, RPA@PBE, G W @PBE, 
and PBE-based rcnormalizcd second-order perturbation 
theor y 62 ' 63 (rPT2) using increasingly larger NAO basis 
sets (Tier 1-4). 

The mean absolute errors of theoretical bond lengths, 
binding energies, and vibrational frequencies as com- 
pared to experiment arc reported in Fig. 1121 The 
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FIG. 12. (Color online) Mean absolute error (MAE) of bond 
lengths (upper panel), binding energies (central panel), and 
vibrational frequencies (lower panel) of LiH, LiF, HF, CO, 
H2, and N2 evaluated at different levels of theory. The esti- 
mated zero-point motion correction has been subtracted from 
the experimental binding energies (reported from Ref. |64T ). 
Calculations were done with a Tier 3 basis set (the largest 
basis set available) for H2, LiH, and HF, whereas a Tier 4 
basis was used for N2 and CO. Numerical values are reported 
in Tables Bill HH and [W] in the Appendix. 



corresponding numerical values are reported in Ap- 
pendix [B] Since our calculations are performed in the 
Born-Oppenheimer approximation with clamped nuclei, 
we compared our calculations with zero-point motion 
corrected experimental binding energies from Ref. |64|. 
PBE, HF, and several perturbative approaches based 
on MBPT, namely G W , EX+cRPA, and rPT2 are in- 
cluded for comparison. 

Self-consistency provides better bond lengths than per- 
turbative GqWq calculations. However, the accuracy 
achieved by sc-GW for the bond lengths is still compara- 
ble to perturbative RPA and not as good as rPT2@PBE, 
which includes higher order exchange and correlation di- 
agrams. 

The binding energies obtained from Go Wo based on HF 
and PBE, reported in Table ITT1 are systematically over- 
estimated. Self-consistent GW over-corrects this trend 
and yields binding energies that slightly underestimate 
experiment. RPA@HF and sc-GW give a very similar 
description of the binding energy, the deviation between 
the two methods being approximately 10 — 20 mcV. This 
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similarity is expected for two reasons: First, in diatomic 
molecules screening is small, thus the sc-GW Green func- 
tion resembles the Hartree-Fock one (since in absence 
of polarization W reduces to the bare Coulomb inter- 
action). Second, the RPA total energy is a variational 
functional of the Green function, and therefore RPA to- 
tal energies are close to sc-GW ones, if the input Green 
function is close enough to the sc-GW Green function^. 
Larger discrepancies between perturbative RPA and sc- 
GW are to be expected for the structural properties of 
systems for which the sc-GW density is substantially 
different as compared to HF or semi-local DFT. Exam- 
ple of these material are molecular interfaces and charge 
transfer compounds, where the ground-state density (and 
the charge transfer) depends strongly on the level align- 
ment between the individual components of the system. 
This will be addressed in future works. Also, for bind- 
ing energies and bond lengths, sc-GW is outperformed 
by rPT2@PBE, which illustrates the importance of in- 
cluding exchange and correlation diagrams beyond the 
GW approximation for a systematic improvement of the 
ground-state properties of finite systems. 

For vibrational frequencies, the dependence on the 
starting point is larger than for binding energies or bond 
lengths. In this case, the best agreement with experi- 
ment is achieved with the PBE functional, whereas for 
HF the errors are substantially larger. Similarly, PBE- 
based Go Wo and EX+cRPA deviate less from experiment 
than HF-based schemes. For instance, the mean absolute 
error of EX+cRPA@HF is approximately a factor of two 
larger that EX+cRPA@PBE, and the same is observed 
for Go Wo. Consequently, sc-GW gives smaller MAEs 
compared to HF-based schemes, but does not improve 
over PBE-bascd perturbative schemes. 



C. Density and dipole moments at self-consistency 

In perturbative approaches, such as Go Wo, EX+cRPA, 
and rPT2, the electron density of a system is defined 
by the eigenstates of the unperturbed reference Hamil- 
tonian - although in principle perturbative corrections 
to the eigenstates of the unpertubed Hamiltonian could 
be calculated. This constitutes a major drawback, as it 
is in part responsible for the well-known starting point 
dependence of perturbative schemes. Self-consistency, on 
the other hand, permits us to incorporate exact exchange 
and dynamical correlation effects into the electron den- 
sity. To illustrate this aspect in sc-GW, we discuss in the 
following the effects of self-consistency on the density and 
the dipole moment of diatomic molecules. 

Figure [TBI shows the density difference between sc-GW 
and HF for the hydrogen fluoride dimer. The density 
difference illustrates the effects of GW correlation on the 
electron density. sc-GW enhances the angular distribu- 
tion of the electron density exhibiting more pronounced 
angular correlation. Moreover, density is shifted from 
the two lobes perpendicular to the molecular axis to the 




FIG. 13. (Color online) Difference between the sc-GW and 
Hartree-Fock densities for the hydrogen fluoride dimer at its 
experimental equilibrium geometry (d = 0.917 A). At self- 
consistency, electron density is shifted from negative (dark) 
to positive regions (light). Units are A -3 and the calculation 
were performed using a Tier 3 basis set. 



TABLE I. Comparison between experimental and theoret- 
ical dipole moments evaluated from sc-GW, PBE, and HF at 
their corresponding equilibrium bond lengths. All values are 
in Debye. 





LiH 


HF 


LiF 


CO 


MAE 


Exp. 


5.88 


1.82 


6.28 


0.11 


- 


sc-GW 


5.90 


1.85 


6.48 


0.07 


0.07 


PBE 


5.63 


1.77 


6.12 


0.20 


0.14 


PBEO 


5.77 


1.81 


6.20 


0.09 


0.06 


HF 


6.04 


1.89 


6.46 


-0.13 


0.17 



bond region, leading to a reduction of the dipole mo- 
ment as compared to HF, which is in better agreement 
with experiment (see Table UJ. 

The dipole moment provides a systematic way to quan- 
tify the quality of the electron density of a system, as it 
is directly comparable with experimental data. We re- 
port in Table UJ the dipole moment of LiH, LiF, HF and 
CO evaluated from sc-GW, PBE, and HF. sc-GW dipole 
moments are in good agreement with experiments and 
reduce the deviation from experiment by approximately 
a factor of two compared to HF and PBE, that tend to 
under and overestimate, respectively. The quality of the 
dipole moment for the small set of molecules presented 
here, indicates that sc-GW is a promising method for the 
description of charge-trasfer compounds, such as molec- 
ular interfaces and hetero-structures. 



VI. CONCLUSIONS 

We have presented an all-electron implementation of 
the fully self-consistent GW method based on a nu- 
meric atom-centered orbital basis in the FHI-aims code^l. 
Self-consistent GW is based on the iterative solution of 
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Hedin's equations with the GW self-energy and polar- 
izability, and therefore it is significantly different from 
partially self-consistent approaches based on perturba- 
tion theory, such as the quasi-particle self-consistent GW 



seller 



-,53.54 



and self-consistency in the eigenvalue: 



,10,66 



In our implementation, the two-particle operators are 
treated efficiently by means of the resolution of identity 
technique^. The introduction of an auxiliary basis for 
the representation of the frequency and time dependence 
of dynamic operators, facilitates an accurate evaluation 
of Fourier integrals that require just few tens of grid 
points. These ingredients allow us to reformulate Hedin's 
equations in a matrix form, that can be solved with stan- 
dard linear algebra packages. 

We presented an assessment of the spectral proper- 
ties of five molecules of interest for organic photo- voltaic 
applications: thiophene, benzothiazole, 1,2,5-thiadiazole, 
naphthalene, and tetrathiafulvalene. For these systems, 
the quasi-particle energies extracted from the sc-GW 
spectral function are found to be in good agreement with 
experimental photo-emission data for all valence states. 
The sc-GW excitation spectrum systematically improves 
over perturbative GqWq based on semi-local DFT and 
Hartree-Fock. This is interpreted as a consequence of 
the mitigation of over- and under-screening errors char- 
acteristic of G Wo@PBE and G W @HF, respectively. 
GqWq based on PBEO, on the other hand, provides re- 
sults in slightly better agreement with experimental data 
than sc-GW. Thus, the PBEO functional appears to be 
close to an optimal starting point for perturbative calcu- 
lations. 

Self-consistent GW total energies based on the 
Galitskii-Migdal formula permit an assessment of 
ground-state and structural properties of molecules. For 
a small set of diatomic molecules we evaluated bind- 
ing energies, bond lengths and vibrational frequencies. 
The bond lengths improve at self-consistency, but still 
have an accuracy comparable to other perturbative meth- 
ods such as exact-exchange with correlation from the 
random-phase approximation. Binding energies are typi- 
cally underestimated compared to experimental reference 
data and do not substantially improve over perturbative 
GqWq calculations. Our results indicate and quantify the 
importance of including vertex corrections - or alterna- 
tively, higher order correlation and exchange diagrams - 
in order to achieve an accurate description of the struc- 
tural properties of molecules. 

Finally, the dipole moments of a set of hctero-atomic 
dimers were studied to investigate the accuracy of the sc- 
GW density. Compared to Hartree-Fock and PBE, the 
sc-GW dipole moments are found in better agreement 
with experiment. These results indicate that sc-GW is a 
promising approach for the description of charge-transfer 
compounds and hetero-j unctions, where the relative or- 
dering of the frontier orbitals influences the charge trans- 
fer at the interface. 
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Appendix A: Rewriting the Galitskii-Migdal formula 

In Hartree atomic units (% = e = m = 1), the elec- 
tronic contribution of the Galitskii-Migdal total energy 



E, 



GM 



i >^ / dx dt lim 
_ J 



9 , ' 

«TTT + ho 

at 



G a (rt,r't' 



(Al) 
Here, ho is the single-particle term of the many-body 
Hamiltonian, i.e., the sum of the kinetic energy operator 
and the external potential due to the nuclei. Introduc- 
ing the equation of motion for the interacting Green's 
function (see e.g. Ref. [671 ) 



- [dr"dt"X a (rt,r"t")G a (r"t",r't') = 

= 6(t - r')5(t - t') , (A2) 



Eq. IA1I can be simplified by eliminating the partial 
derivative with respect to time, obtaining: 



E GM = -i ]T] I dv dr'dt dt' 



lim [(-V 2 r + 2 Vcxt (r) + m (r)) 8(r - r')S{t - t') 



r — >r 
t'-yt+ 



£ CT (r£,rY)]G CT (rY,r , i') 



(A3) 



Making use of the matrix representation of the Green 
function G CT (r,r',T) = ^\j ¥i( r )Gij ( T )<Pj( T ') ~ with 
G = s~ 1 G a s~ 1 - the first three terms in Eq. IA3I can 
be rewritten as: 

-i^2 dr }}™ r [~ V r + 2u c*t(r) + va(r)] x 
x^ i (r)G;.(7 = 0> i (r') = 

ij 

-*EE^/( r = (r )[ 2 ^ + 2 ^ t +^] • ( A4 ) 



a ij 
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For time independent Hamiltonians, the Green func- 
tion depends only on time differences r = t — t'. In 
the last step of Eq. IA41 we defined the matrix rep- 
resentation of the kinetic energy operator as Uj = 

J dvtpi (r) 

for 



2 



y5j(r), and use a similar representation 

and v^. Finally, the last term in Eq. IA3I 
can be rearranged by using the Fourier transform 
of the Green function and the self-energy G(t,t') = 
I- °° l^e -1 ^* - * ^G(uj), and substituting the matrix rep- 
resentation of G: 

-iV fdrdr"dtdt" lim S <T (ri,r"t")G" T (r"t", r'i') = 

J r'->r 

-EE/^^h^h^ • ( A5 ) 



a ij 



Summing Eqs. lA4l and lA5l and separating the self-energy 
in its exchange and correlation components Ejj(u;) = 



T,fj +Sy(a;), one finally arrives at the expression reported 
in Eq. [301 We refer to Ref. [68|, for a discussion on the 
evaluation of the Galitskii-Migdal formula on the imagi- 
nary frequency axis. 



Appendix B: Binding energies, bond lengths, and 
vibrational frequencies of diatomic molecules 

In Tables HI1 and IIIII we report counterpoise corrected-^ 
binding energies and bond lengths for H 2 , LiH, LiF, HF, 
N2, and CO. The corresponding vibrational frequency are 
reported in Table HVl The sc-GW results are compared 
with experimental value s 64 ' 65 and several perturbativc 
approaches based on MBPT, namely G W , EX+cRPA, 
and rPT2@PBE£2. HF and PBE are included for com- 
parison. 
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TABLE II. sc-GW , and perturbative Go Wo binding energies (evaluated from the Galitskii-Migdal formula) of diatomic 
molecules compared to (zero point motion corrected) experimental reference data taken from Ref. [6J. We report pertur- 
bative RPA, HF, PBE, and renormalized second-order perturbation theory (rPT2) for comparison. Calculations were done 
with a Tier 3 basis set (the largest basis set available) for H2, LiH, and HF, whereas a Tier 4 basis was used for N2 and CO. 
The mean absolute errors (MAE) are reported in panel (b) of Fig. 1121 All values are in eV. 





H 2 


LiH 


HF 
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PBE 
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TABLE III. sc-GW and perturbative Go Wo bond lengths of diatomic molecules compared to experimental reference data taken 
from Ref. [63. RPA, HF, PBE, and renormalized second-order perturbation theory (rPT2) are included for comparison. Mean 
absolute errors (MAE) are reported in panel (a) of Fig. 1121 All values are in A. 





H 2 
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N 2 
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TABLE IV. sc-GW , and per turbative GoWq vibrational frequencies of diatomic molecules compared to experimental refer- 
ence data taken from Ref. |6g. RPA, HF, PBE, and renormalized second-order perturbation theory (rPT2) are included for 
comparison. The mean absolute errors (MAE) are reported in panel (c) of Fig. 1121 All values are in cm -1 . 
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